###############################################################################
# Function that return slopes in a data table (coef and se)
###############################################################################

slopes_dt = function(f,names_var) {
  nn = rownames(f$coeftable)
  slopes_pos = grep("ddr:",nn)
  dt_plot = data.table(f$coeftable[slopes_pos,])
  nn_pos = substr(nn[slopes_pos],17,30)
  nn_pos = as.numeric(nn_pos)
  dt_plot=dt_plot[,1:2]
  dt_plot[,fips:=nn_pos]
  setnames(dt_plot,c(paste0(c("coef","se"),names_var),"fips"))
  return(dt_plot)
}

###############################################################################
# Function that return slopes in a data table (coef and se) when slopes are 
# per periods
###############################################################################

slopes_dt_period = function(f,np) {
  nn = rownames(f$coeftable)
  
  # First period
  name_now=paste0("ddr_p",1,":")
  slopes_pos = grep(name_now,nn)
  dt_plot = data.table(f$coeftable[slopes_pos,])
  nn_pos = substr(nn[slopes_pos],20,30)
  nn_pos = as.numeric(nn_pos)
  dt_plot=dt_plot[,1:2]
  dt_plot[,fips:=nn_pos]
  setnames(dt_plot,c("coef","se","fips"))
  dt_plot[,period:=1]
  dt_tot = copy(dt_plot)

  # For the following periods
  for (k in seq(2,np)) {
    name_now=paste0(":ddr_p",k)
    slopes_pos = grep(name_now,nn)
    dt_plot = data.table(f$coeftable[slopes_pos,])
    nn_pos = substr(nn[slopes_pos],13,17)
    nn_pos1 = as.numeric(nn_pos[!is.na(as.numeric(nn_pos))])
    nn_pos2 = as.numeric(substr(nn_pos[grep(":",nn_pos)],1,4))
    nn_tot =c(nn_pos2,nn_pos1)
    dt_plot=dt_plot[,1:2]
    dt_plot[,fips:=nn_tot]
    setnames(dt_plot,c("coef","se","fips"))
    dt_plot[,period:=k]
    dt_tot = rbind(dt_tot,dt_plot)
    
  }
  return(dt_tot)
}


###############################################################################
# Function that generates and saves table
###############################################################################

sum_table = function(data,name_vars,save_name="",weight_var) {
  
  dt=copy(data)
  dt[,.w:=.SD,.SDcols=weight_var]
  
  # Percentiles, mean and variances
  s=as.matrix(dt[,lapply(.SD,wtd.quantile,weights=w,probs=c(0.1,0.25,0.5,0.75,0.9),na.rm=TRUE),
                       .SDcols=name_vars])
  s=rbind(s,as.matrix(dt[,lapply(.SD,wtd.mean,weights=w,na.rm=TRUE),.SDcols=name_vars]))
  s=rbind(s,as.matrix(dt[,lapply(.SD,wtd.var,weights=w,na.rm=TRUE),.SDcols=name_vars]))

  rownames(s)=c(paste0("Percentile ",c(10,25,50,75,90)),"Mean","Variance")
  l=dim(s)[2]
  dig=rep(3,l+1)
  if (nchar(save_name)>1) {
    print(xtable(s,type="latex",digits=dig,display=c("s",rep("f",l))),hline.after = NULL,
          file=paste0(paste0(save_results,save_name),".tex"),include.rownames = TRUE,include.colnames=FALSE,
          sanitize.text.function = function(x){x},only.contents = TRUE)
  }
  
  return(s)
}


###############################################################################
# Function for maps
# (needs dt_map in the environment)
###############################################################################

plot_map = function(data,bmin=0,bmax=0,trunc=0,save_name="",llabel="Temperature effects     ",palette="Spectral",nc=11,ingrey=FALSE,rev=0) {
  
  if (ingrey==TRUE & rev==0) {
    #colors_now = gray.colors(n=nc,rev=FALSE)
    colors_now = gray.colors(n=nc,start=0.1,end=0.9,rev=FALSE)
    #colors_now = gray.colors(n=nc,start=0.15,end=0.85,rev=FALSE)
    color_bound = "white"
  } 
  if (ingrey==TRUE & rev==1) {
    #colors_now = gray.colors(n=nc,rev=FALSE)
    colors_now = gray.colors(n=nc,start=0.1,end=0.9,rev=FALSE)
    #colors_now = gray.colors(n=nc,start=0.15,end=0.85,rev=TRUE)
    color_bound = "white"
  } 
  if (ingrey==FALSE) {
    color_bound = "gray90"
    myPalette <- colorRampPalette((brewer.pal(nc, palette)))
    myPaletteR <- function(x) rev(myPalette(x))
    if (rev==0) {
      colors_now = myPalette(100)
    } else {
      colors_now = myPaletteR(100)
    }
  }
  
  dt_groups=copy(data)
  setnames(dt_groups,c("fips","GFE"))
  
  # Truncate values to bmin and bmax
  if (trunc!=0) {
    dt_groups[GFE< bmin,GFE:=bmin]
    dt_groups[GFE> bmax,GFE:=bmax]
  }  
  if (trunc==0) {
    bmax = max(abs(min(dt_groups$GFE)),abs(max(dt_groups$GFE)))
    bmin=-bmax
  }
  
  # # Add the fips with missing values
  # lf = unique(dt_map[is.na(fips)==FALSE,fips])
  # dim(dt_groups)
  # dt_groups = merge(dt_groups,data.table(fips=lf),by="fips",all=TRUE)
  # dim(dt_groups)
  # summary(dt_groups)
  
  # # Add variable of missing mg
  # dt_groups[,miss:=is.na(GFE)]
  
  # Merge data with map
  dt_plot = merge(dt_map,dt_groups,by="fips",all.x=T)
  summary(dt_plot)
  
  # Keep only counties with values
  length(unique(dt_plot[is.na(GFE)==FALSE,fips]))
  dt_plot = dt_plot[is.na(GFE)==FALSE]
  length(unique(dt_plot[is.na(GFE)==FALSE,fips]))
  
  # Plot
  g=ggplot(data=dt_plot[long>-100],aes(x=long,y=lat,group=group, fill=GFE)) +
    geom_polygon(color = color_bound, linewidth = 0.01) +
    coord_map(projection = "albers", lat0 = 45, lat1 = 55) +
    labs(fill=llabel) +
    theme(legend.position="bottom",
          axis.line=element_blank(),
          axis.text=element_blank(),
          axis.ticks=element_blank(),
          axis.title=element_blank(),
          panel.background=element_blank(),
          panel.border=element_blank(),
          panel.grid=element_blank(),
          legend.key.width=unit(1.2,"cm"),
          legend.text = element_text(size=18),
          legend.title = element_text(size=18))+
    guides(fill = guide_colorbar(barwidth = unit(1.8, "inches"),  # bar width
                                 barheight = unit(0.2, "inches")))  
  
  sc = scale_fill_gradientn(colours = colors_now, limits=c(bmin,bmax))
  g1= g+sc
  print(g1)
  
  if (nchar(save_name)>1) {
    
    #ggsave(paste0(save_res,"/",save_name,".pdf"))
    ggsave(paste0(save_res,"/",save_name,".pdf"),width=7, height=4.95,dpi=300)
  }
}
